
    ##########################################################################################
    ### R code for Figure 6 of Supplemental Material
    ### 
    ##########################################################################################
 
    ### Note: This code requires "Summary_Outputs.csv" and "Bill.Param.csv", which will be produced after running updated "Figure3_R_code.R". 
    
    DW_NOM_Rates <- read.csv("Supplemental_Figure6_Data1.csv") # DW-NOMINATE Classification Success Rates for Foreign Policy Votes (Sorted from Voteview.com)
    Ideals <- read.csv("Summary_Outputs.csv")
    Bill.Param <- read.csv("Bill.Param.csv")

    T <- 33  # number of congresses
    year.odd <- seq(from=1945, to=2010, by=2)


    ##############################################################
    ### Get the classification rates
    class_rates <- matrix(NA, nrow=T, ncol=2)

    for (t in 1:T){

    # merge Vote_Data with Ideal Point Estimates
    Vote_Data <- read.table(paste("All_FP_votes_s", t+78, "th.txt", sep=""), header=T)
    Ideal <- subset(Ideals, time==t, select=c(ICPSR_ID, time, congress, IRT_mean))
    Data <- merge(Ideal, Vote_Data, by="ICPSR_ID", all=F)
    x.mean <- Data$IRT_mean
    N <- length(x.mean)
    Y <- Data[,-c(1:8)]

    Bill.p <- subset(Bill.Param, time==t, select=c(time, congress, vote.no, mean_alpha, mean_beta))
    alpha <- Bill.p$mean_alpha
    beta  <- Bill.p$mean_beta
    M <- length(alpha)

    # (1) Classification Rate
    P <- matrix(NA, N, M)

    for (i in 1:N){  # Loop over legislators
    for (j in 1:M){
            P[i,j] <- plogis(x.mean[i]*beta[j] - alpha[j])
    }
    }

    ### Generate Y.pred

    Y.pred <- matrix(NA, N, M)
    for (i in 1:N){
    for (j in 1:M){
    if (P[i,j] >= .5) Y.pred[i,j] <- 1
    else Y.pred[i,j] <- 0
    }}

    ### Clean up data
    Y1 <- as.matrix(Y)
    Y1[is.na(Y1)] <- 9 # To change "NA" to "9". 

    # 1. For all legislators
    Count <- matrix(NA, N, M)
    for (i in 1:N){
    for (j in 1:M){
    if (Y1[i,j] == 9) {Count[i,j] <- 9}
    else {if (Y1[i,j] == Y.pred[i,j]) {Count[i,j] <- 1}
          else Count[i,j] <- 0}
    }}
    Count1 <- Count
    Count1[Count1==9] <- NA

    # to calculate the total number of NA
    NA.count <- matrix(0, N, M)
    NA.count[is.na(Count1)] <- 1

    # to calculate the total number of 1's (=right classification)
    one.count <- matrix(0, N, M)
    one.count[Count==1] <- 1

    # to calculate the total number of 0's (=wrong classification)
    zero.count <- matrix(0, N, M)
    zero.count[Count==0] <- 1

    # Classification rate for all
    Class.All <- 100*(sum(one.count)/(sum(one.count)+ sum(zero.count)))

    class_rates[t,2] <- Class.All
    class_rates[t,1] <- t+78
    }

    ### Now, Compare IRT and DW-NOM
      
    postscript("Supplemental_Figure6.eps", height=6, width=8, horizontal=F)
    par(mfrow=c(1, 1))
    plot(1:T, class_rates[,2], type='o', lty=1, col='blue', xaxt="n", xlab="", ylab="Overall Classification Success Rate (%)", ylim=c(50, 100))
    points(1:T, DW_NOM_Rates[,2]*100, type='o', lty=2, col='red')
    axis(side=1, las=3, at=1:T, labels=year.odd, cex.axis=.8)
    #title(main="The Numbers of Moderate Senators by Party",  font.main=1)
    legend('topleft', legend=c("Bayesian IRT", "DW-NOMINATE"), lty=c(1, 2), col=c('blue', 'red'), cex=.8)
    dev.off()
    


